Retracing cyanobacteria blooms in the Baltic Sea

In late summer, massive blooms and surface scums of cyanobacteria emerge regularly in the Baltic Sea. The bacteria can produce toxins and add bioavailable nitrogen fixed from atmospheric nitrogen to an already over-fertilized system. This counteracts management efforts targeted at improving water quality. Despite their critical role, the controls on cyanobacteria blooms are not comprehensively understood yet. This limits the usability of models-based bloom forecasts and projections into our warming future. Here we add to the discussion by combining, for the first time, satellite estimates of cyanobacteria blooms with output of a high-resolution general ocean circulation model and in-situ nutrient observations. We retrace bloom origins and conditions by calculating the trajectories of respective water parcels backwards in time. In an attempt to identify drivers of bloom development, we find that blooms originate and manifest themselves predominantly offshore where conditions are more nutrient-depleted compared to more coastal environments.

Reports of cyanobacteria blooms all over the world have been increasing over the years 1 . This also applies to the Baltic Sea where actual trends in late-summer bloom occurrences (predominantly Nodularia spumingena) are, however, discussed controversially [2][3][4][5][6][7] . Related to such cyanobacteria blooms, toxins might enter into the food web 8,9 and bioavailable nitrogen (bioavailable N) is added to the system by their ability to fix atmospheric nitrogen 10 . This input by nitrogen fixation is substantial (up to 40% compared to the overall nitrogen loads from rivers and atmospheric deposition 11,12 ) and fertilises additional primary production and associated export of organic matter to depth. The latter, in turn, is associated with increased oxygen consumption at depth. The related oxygen depletion has been shown to release phosphorus from anoxic sediments 13,14 . This can alleviate phosphorus limitation of autotrophic growth at the sun-lit surface as mixing and upwelling closes the positive feedback loop which may accelerate deoxygenation.
In summary, there is concern because cyanobacteria add nutrients to an already over-fertilized ecosystem thereby degrading the water quality 15,16 and there is the apprehension that this process may accelerate as we move further into the Anthropocene 17 . Thus, efforts targeted at a comprehensive understanding are increasing as the societal relevance becomes more apparent 8,[18][19][20][21][22] . Even so, several mechanisms behind the formation of a cyanobacteria bloom are yet to be deciphered 23,24 and reliable quantitative tools for projections are still under development 25 . We argue that the former must be accomplished in order to reliably assess the role of potentially antagonistic anthropogenically-induced eutrophication 14,26,27 and climate change 28 .
To date, numerous abiotic and biotic factors have been suggested to promote cyanobacteria blooms [29][30][31][32][33] . As concerns their implementations in numerical models the underlying parametrizations are in large parts based on empirical field correlations in a process where modellers adjust poorly constrained model parameters (e.g. the saturation coefficients for light and nutrients) such that model output aligns with (typically noisy) observational data. Among the positive correlations exploited so far (for e.g. Nodularia spumigena) are: (1) weather conditions and photosynthetically active radiation, (2) ambient temperatures exceeding 16 • C and (3) elevated phosphate concentrations concurring with relative nitrate depletion 29,30 . Negative correlations include salinities exceeding 10 PSU 32 and winds in excess of a moderate breeze. In this context we face an apparent paradox where relatively slow-growing diazotrophs compete for bioavailable phosphorus with faster-growing phythoplankton -while they, ultimately stock up the system with nitrogen that is essential to their competitors which can not utilise atmospheric dinitrogen [34][35][36][37] . In numerical models the extinction of diazotrophs (unless downstream of denitrification zones 38 ) is prevented by diverse approaches which include factoring in diazotrophic capability to utilize dissolved organic phosphorus (DOP) effectively 39 , to utilize relatively low P-concentrations 40 , and selective grazing where cyanobacteria experiences less grazing pressure than other phytoplankton [41][42][43] .
The present study adds to the ongoing discussion by merging satellite data, high-resolution ocean circulation modelling and in-situ observations to dissect the bloom-triggering factors. The aim is to overcome problems www.nature.com/scientificreports/ associated with the (limited) availability of observational data and to include not only the bloom peak but also the bloom formation phase into our analyses. Blooms are highly variable in space and time which translates into a relative sparseness of observational data despite substantial observational effort. Consequently, model evaluations of simulated diazotrophy in the Baltic Sea have mostly been limited to long-term averages or interannual variations of cumulative proxies indicative of standing stocks 43,44 . Our focus is on environmental conditions preceding bloom manifestations (as picked up by satellite) rather than on conditions during blooms (which were already investigated in the past 45,46 ). Essentially, this study compares trajectories of water parcels that develop blooms with those that don't in the Baltic Proper.

Results
Tracing the bloom origins. Figure 1a shows the initialisation for the bloom event identified from space on July 11th, 2010 52 . Randomly seeded drifters (i.e., randomly distributed virtual floats) were attributed to blooming and non-blooming water parcels based on satellite information. Drifters that were seeded at locations of inconclusive satellite information were discarded (further details on this procedure are provided in the Method section). Questions we set out to answer are: What are the differences among bloom-forming water parcels and those that develop no bloom? Do they share a common history of environmental conditions or can we dissect differences hinting at controls of diazotrophy?
History of water parcels. The history of water parcels is retraced using the high-resolution general ocean circulation model MOMBA 47,48 . Starting from initial (random) positions trajectories are calculated for 60 days backwards in time using archived MOMBA output of surface currents. We implicitly assume that all blooming is preconditioned within 2 months. Figure 1b exemplarily shows the retraced origin of each water parcel 30 days prior to the blooming event on July 11th, 2010. Respective animations for all events considered here (one per year throughout 2010 to 2015 as shown in Fig. 2) are archived at http:// www. baltic-ocean. org/ categ ory/ cyano. Figure 3 summarises the abiotic conditions to which water parcels are exposed throughout the 60 days leading up to a blooming event: Contrary to our expectations we can not identify robust links between sea surface temperatures, photosynthetically active radiation and surface mixed layer depth, on the one hand, and blooming on the other hand. The respective differences among parcel trajectories leading to blooms, and trajectories which do not, are very small and appear insignificant. For example, sea surface temperatures associated with blooms are on average merely 0.25 • colder and do not occur consistently in all of the considered years. Further, photosynthetically active radiation and surface mixed layer do also fail to exhibit consistent differences throughout the 60 day periods.
The surprise, regarding these very small temperatures and mixed layer depths differences, is that cyanobacteria have been shown to grow more optimal in warmer temperatures 25 and, further, that the high energetic requirements of utilising the exceptionally chemically stable atmospheric dinitrogen molecules suggest a preference for higher radiation levels as given by, e.g., shallower surface mixed layers pervaded by solar radiation (and higher radiation levels). While abiotic environmental conditions such as radiation and temperature that are known to affect physiological processes in cyanobacteria cells 25 lack all explanatory power to distinguish between blooming and non-blooming water parcels, the distance to the coast is apparently strongly linked to blooming probability: Fig. 3d shows that water parcels that result in blooming separate from the coast and travel offshore. The separation starts weeks prior to blooming which allows enough time for differing environmental conditions to manifest themselves in blooming dynamics. More specifically we find that 90% of blooming water  Typical nutrient characteristics. The biochemical properties of near coastal waters are known to differ from the waters further off-shore 50 . Typically they are more enriched in nutrients that are essential for algae growth because coastal upwelling taps into deep nutrient replete waters and thereby fertilises the sun-lit surface ocean. Figure 4 confirms that offshore waters that are more prone to develop cyanobacteria blooms are characterised by more oligotrophic, nutrient-depleted conditions than the more coastal waters that bloom less: all available (irregular sampled) ICES surface nutrient observations from 1990-2017 are binned into climatological months. Based on our results from Fig. 3 we define coastal distances of more than 60 nautical miles as "off-shore" and distances closer than 40 nautical miles as near "coastal" and find, amidst great variability, that offshore waters are generally more nutrient-depleted than waters closer to the coast. This applies both to nitrate and phosphate surface concentrations except for very early in May.

Discussion
We set out to find links between environmental conditions and the occurrences of cyanobacteria blooms. Foregoing studies focussed on conditions during the peak of the blooms 45 , or on the frequency of bloom occurrences over the years 46 . In contrast, our approach of combining high-resolution ocean circulation modelling, satellite data, in-situ observations and a back tracing algorithm made it possible to compare the history and trajectories of water parcels comprising a later bloom with the history of "non-blooming" parcels. To this end we add to the discussion on controls of diazotrophy by providing an analysis that is inherently free of seasonal biases and limelights conditions during the growth phase of cyanobacteria in comparison to non-bloom forming waters. Our results differ from earlier findings that stress the importance of photosynthetically available (solar) radiation and surface mixed layer depth (or wind) for bloom formation. We could also not confirm the hypothesis that blooms may be triggered by strong wind-induced mixing events 51 -even though our approach was designed to confirm and quantify just these links. Moreover, we found no indication of an influence either of insolation or of temperature other than that blooming occurs generally in summer.
The only apparent difference between bloom forming waters and other water masses in our analysis is the distance to the coast: trajectories of water parcels that develop a bloom are typically situated further offshore than those that do not. Because coastal and offshore waters have different biogeochemical characteristics 50 this www.nature.com/scientificreports/ may provide a clue as for what may trigger the blooms. Specifically, we find that average nutrient concentrations of phosphate and nitrate are lower offshore. This applies also to P * (defined as DIP in excess to DIN/16). Our results apparently contradict a common assumption that P * is a major driver of cyanobacteria blooms in the Baltic and we speculate that ecological complexity obscures simple links between physiological traits and manifestations of blooms in field data.
Caveats apply. Among them are uncertainties associated with the satellite retrieval algorithm and the realism of the ocean-circulation model used to back-trace the conditions that led to blooming events.  (Fig. 2). The onset of these blooms is generally relatively sudden with only smaller stray bloom patches in the area some days in advance (cf. respective visualizations archived at http:// www. baltic-ocean. org/ categ ory/ cyano).

Methods
After identifying the blooms in satellite observations backward trajectories up to 2 months prior to bloom occurrence were calculated. Initial positions of 4000 drifters were randomly seeded within a region bounded by 13 • E to 24 • E and 54 • N to 60 • N. These drifters were then attributed to blooming and non-blooming patches. Patches where the satellite data were inconclusive were discarded. Conditions prevailing within these parcels during their 2 month voyage through the Baltic were analysed and compared for bloom forming and non-bloom Nutrient observations. In-situ observations of surface nutrients were provided by the International Council for the Exploration of the Sea (ICES) with a precision of 10 −2 for phosphate and 10 −1 for nitrate. Despite the extraordinarily dense data coverage (when compared to typical open ocean conditions e.g. in the Atlantic), the data were still too sparse for robust conclusions based on measurements that were taken in the water parcels that were actually tracked. We mitigated this problem by binning data into: (1) observations during summer seasons 1990-2017 and (2) into near coastal and offshore data. Hereby coastal was defined as waters less than 40 nautical miles off the coast while offshore values refer to all measurements more than 60 nautical miles off the coast. This choice is motivated by Fig. 3d. Overall our approach results in 8146 coastal and 3819 offshore observations   N). The model features a horizontal resolution of ≈1.9 km (corresponding to one nautical mile) and is eddy-rich in that it starts to resolve the relevant spatial scales 54 . The vertical discretization comprises 47 levels.
There are no open boundaries and the model domain is surrounded by solid walls. We use the K-profile parameterization (KPP 55 ) with parameters identical to those applied in other eddy-permitting configurations 56,57 . The difference in MOMBA settings relative to 47 is the atmospheric boundary condition which changed to reanalysis data provided by the Swedish Meteorological and Hydrological Institute based on the atmospheric model RCA4 58 .
Backward trajectories. In order to identify the conditions that lead to cyanobacteria blooms (in relation to conditions that do not) we assume that the flow field consists of (Lagrangian) water parcels, each associated with a volume of surface water. 4000 water parcels for analysis are randomly chosen (or seeded) within the central Baltic Sea (ranging from13 • E to 24 • E and 54 • N to 60 • N) during the bloom peaks in 5 consecutive years, identified as described in Sect. 2.1.1. The parcels (or virtual drifters) are divided into two groups -parcels that are located within identified cyanobacteria bloom patches and those who are outside blooming patches. (A third group related to patches where the satellite classification algorithm is inconclusive is discarded). An example of parcel allocation for a bloom in 2010 is depicted in Fig. 2. The trajectory of each water parcel is calculated backward in time from the simulated Eulerian surface velocity fields as calculated by the ocean circulation model MOMBA based on 3-hourly simulated surface velocity model output. For all calculations, the velocities, that are defined on the MOMBA spatial model grid are linearly interpolated in space and time onto the (moving, backtraced) positions of water parcels. Back-tracing is calculated up to 60 days prior to respective blooming using a generic Euler forward (or rather backward) algorithm. Note, that a similar approach to calculate backward trajectories has been used and tested earlier for other purposes and other regions 59,60 .
In a second step, the back-traced trajectories are linked to their respective simulated abiotic water mass properties (such as surface water temperature, mixed layer depth and photosynthetically active radiation (RAR) within the mixed layer). PAR is calculated from the incoming solar radiation averaged over the mixed layer while assuming that PAR provides 42% to the incoming solar radiation and assuming a constant light attenuation of 0.2 m −161 . As a caveat we report here that shading effects of phytoplankton blooms are not included in this calculation. This effect has been suggested to be important e.g. in lakes whereby cyanobacteria evade shading-effects of ordinary phytoplankton by increasing their buoyancy 62 . In the Baltic Sea, cyanobacteria surface-scums are presumably the major plankton-induced shading events. The latter renders an inclusion of shading effects problematic because causal relationships may be obscured if the presence of cyanobacteria is linked to environmental conditions caused by themselves.

Data availability
Movies containing the satellite data and the backtraced drifters are archived in form of movies at http:// www. baltic-ocean. org/ categ ory/ cyano. The tracking data are available at https:// doi. org/ 10. 5281/ zenodo. 66415 49. The data are distributed under the Creative Commons Attribution 4.0 License.